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ABSTRACT 

In this letter, we describe results of new high-resolution axisymmetric relativistic 
MHD simulations of Pulsar Wind Nebulae. The simulations reveal strong breakdown 
of the equatorial symmetry and highly variable structure of the pulsar wind termina- 
tion shock. The synthetic synchrotron maps, constructed using a new more accurate 
approach, show striking similarity with the well known images of the Crab Nebula 
obtained by Chandra, and the Hubble Space Telescope. In addition to the jet-torus 
structure, these maps reproduce the Crab's famous moving wisps whose speed and 
rate of production agree with the observations. The variability is then analyzed using 
various statistical methods, including the method of structure function and wavelet 
transform. The results point towards the quasi-periodic behaviour with the periods of 
1.5 — 3yr and MHD turbulence on scales below lyr. The full account of this study 
will be presented in a follow up paper. 
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1 INTRODUCTION 

The Crab Nebula is the archetypal compact synchrotron 
nebula, continuously powered by a magnetized, ultrarela- 
tivistic wind from a young, rapidly rotating neutron star. 
In order to match the sub-relativistic expansion of the con- 
fining medium (either interstellar medium (ISM) or super- 
nova remnant (SNR)), the wind must be decelerated. This 
occurs at a strong termination shock, established at a ra- 
dius where the ram pressure of the pulsar wind matches 
the pressure of its surrounding ( Ree s &; Gunn|1974 ; Kennel 
& Coroniti||1984a|). This shock causes heating of the wind 
plasma and presumably particle acceleration. As a result a 
bubble of relativistic particles and magnetic field, which is 
known as a pulsar wind nebula (PWN), is formed between 
the termination shock and the supernova shell (for general 
reviews see |Gaensler & Slane| (2006); |Bucciantini| ([2008a); 



Kirk, Lyubarsky & Petri (2007)). This way, the spin-down 
energy lost by the pulsar through its wind becomes de- 
tectable as non-thermal radiation, mostly synchrotron. 

The steady progress of X-ray astronomy has resulted in 
the discovery of a rather peculiar jet-torus structure in the 
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inner part of the Crab Nebula, not predicted by the classical 
model of plerions. This structure had been seen in earlier 



observations by Brinkmann et al. (1985) and Hester et al. 
( 1995), but the recent observations using the Chandra X-ray 
satellite and the Hubble Space Telescope show the structure 
and its temporal nature in spectacular detail (Weisskopf et 



al.|20 00; Hest er et al.|2002 ). Similar structures are found in 



many other PWNe (|Weisskopf et al.|2000||Slane et al.|2 004 



Romani & Ng 2003; Camilo et al. 2004; Gaensler et al. 2001 
2002| |Gotthelf fc Wang|2000| |Helfand et al.|2001| |Hester et 



al.||1995||Lu et al.|2002| |Pavlov et al"[ 2003 ) . 

The first theoretical models for the formation of both 
the equatorial torus and the polar jets, that combine to form 
the jet-torus, were presented by Bogovalov & Khangoulian 
( 2002p and Lyubarsky (2002). Both these models agree in 
that the torus arises due to the anisotropic distribution of 
the pulsar wind power, enhanced in the equatorial direction, 
and the jets are formed downstream of the wind termination 
shock. In particular, Lyubarsky (2002) proposed that jets 



are produced via the so-called "tooth-paste" effect related 
to the hoop stress of the azimuthal magnetic field carried 
by the wind. These theoretical models were confirmed, in- 
cluding the "tooth-paste" effect, by series of recent numer- 



ical simulations (Komissarov & Lyubarsky 20031 2004 Del 



Zanna et al. 2004 Bogovalov et al. 2005}, performed us- 
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ing newly developed codes for Relativistic MHD (RMHD). 
The RMHD model of PWNe has been recently developed 
to a high degree of sophistication, and has proven itself ca- 
pable of explaining, at least qualitatively, many properties 
of PWNe, including their morphology, polarization, spec- 



tral index maps and gamma-ray emission ( Bucciantini et al 



2005| |Del Zanna et al.|| 2006| |Volpi et al.||2008[ |2009 ) 

Although the global structure of Crab's jet-torus ap- 
pears to be quite persistent over long time-scales, it has 
been known since 1920s that its fine structure shows a short 
time variability both in the vicinity of the termination shock 
and at larger distances ( |Lampland|192l[ [Scargle 1969). The 
most fascinating variable features are the thin synchrotron 
emitting filaments, dubbed by I Scargle (1969) as wisps, who 



concluded that they may move outwards with relativistic 
speeds. However, it is the recent high resolution optical and 
X-ray monitoring observations of the Crab Nebula that have 
enabled us to firmly establish that this is indeed the case 
( Tanvir et al.||1997| [Hester et al.||2002| ). They have also dis- 
covered the jet variability. In particular, the proper motion 
of the jet emission pattern clearly indicates relativistic out- 
flow. Moreover, rather bright and variable features appear 



near the jet base, like the so-called sprite ( Hester et al. 2002 ) 



Scargle (1969) commented that the wisps could be ei- 



ther 1) inhomogeneities advected by the flow in the nebula 
or 2) magnetosonic waves. Later, Gallant & Arons (1994) 



noted that, if ions are part of the pulsar wind, then down- 
stream of the termination shock their Larmor radius can be 
comparable to the separation between the brightest wisps 
of the Crab Nebula, and thus the wisps may also reflect the 
steady kinetic structure of the termination shock as medi- 
ated by the ions. This idea was recently developed further 



by Spitkovsky & Arons (2004) who have shown that the 



cyclotron instability of ions may lead to the excitation of 
magnetosonic waves in the electron-positron component of 
the post termination shock flow and thus explain the ob- 



served variability of the wisps. Begelman (1999) proposed 
an alternative, purely MHD model of the wisp origin. He 
assumed that the flow behind the termination shock con- 
sists of a fast equatorial zone and a slower moving surround- 
ing zone, and that the shear between the zones gives rise to 
the Kelvin-Helmholtz instability which causes the equatorial 
zone to ripple, these moving ripples are identified with the 
dynamic wisps of the Crab Nebula. Recently, Bucciantini & 
|Del Zanna] ( |2006| ) argued that the growth rate of the Kelvin- 
Helmholtz instability is too slow to explain the wisps in this 
way. On the other hand, the previous RMHD simulations 
of PWN did indicate noticeable variability of the flow and 
suggested the possibility of explaining the wisp production 
in the MHD framework (|Komissarov &; Lyubarsky 2004 Bo- 



govalov et aLp005l |Volpi et al.|[2 008 ) 

In order to explore whether the observed variability of 
the Crab Nebula can indeed be explained within the purely 
magnet ohydro dynamic approach we have carried out new 
RMHD simulations with the highest numerical resolution 
so far, and developed new sophisticated method for the con- 
struction of synthetic maps of the synchrotron emission from 
the simulated PWN. In addition, we did not impose the re- 
flective boundary in the equatorial plane. The simulations 
have revealed the highly unsteady non-linear behaviour of 
the termination shock and the highly intricate fibrous struc- 
ture of the PWN emission, including bright moving wisps. 



Here we briefly describe the key results of this study whereas 
the full account will be given elsewhere. 



2 DETAILS OF NUMERICAL SIMULATIONS 

To perform these simulations we used an improved version 



of the RMHD scheme constructed by Komissarov (1999) 



In order to reduce numerical diffusion we apply parabolic 
reconstruction instead of the linear one of the original code. 
Our procedure, in brief, is to calculate the minmod- averaged 
first and second derivatives and use the first three terms of 
the Taylor expansion for spatial reconstruction. 



2.1 Computational grid and initial solution 

We utilise a spherical grid with standard coordinates {r, 
where in the direction we adopt a fixed angular cell size AO. 
In the radial direction we adopt a logarithmic scaling, with 
radial cell spacing growing as An = riAO. As the Courante 
stability condition requires At < An/c, we split the com- 
putational domain into a set of rings such that the outer 
radius of each ring is twice its inner radius, and advanced 
the solution in each jth ring with with its own time step Atj , 
where Atj+i = 2 Atj. In order to ensure that the termina- 
tion shock is always fully inside the computational domain 
and to reduce the computational cost, which is dominated 
by the inner rings, the number of rings is variable and is kept 
to the minimum needed to capture the termination shock. If 
the termination shock moves very close to the inner bound- 
ary an additional inner ring is activated. When it moves 
outside of the inner ring, this ring is deactivated. Moreover, 
in order to reduce the computational cost further more, we 
do not update the cells which are fully inside the pulsar wind 
zone. In order to verify the convergence of numerical solu- 
tions, in the statistical sense, we used four grids with 100, 
200, 400, and 800 cells in the direction. In this letter, we 
only present the results obtained with the highest resolution, 
which exceeds by the factor of 4 the resolution of previous 
simulationsQln the radial direction the computational do- 
main extends up to the distance corresponding to r = lOlyr 
and the simulations are run up to the time corresponding to 
the current age of the Crab Nebula, t ~ 960 yr. 

The initial pulsar wind zone is spherical with the ra- 
dius r w = 0.5 lyr. The wind model is exactly the same as 



in Komissarov & Lyubarsky (2004) with the magnetization 
parameter £ = 0.4, the Lorentz factor W = 10, and purely 
azimuthal magnetic field changing sign at = tt/2 accord- 
ing to the model of dissipation of the alternating component 
of magnetic field of striped wind at the termination shock 
by Lyubarsky (2003b). Initially, the region outside of the 



wind zone is filled with cold radially expanding flow with 
constant density p e and velocity v oc r. The parameters of 
this flow are selected in such a way that its total mass and 
kinetic energy within the radius r e = 1.47 lyr are 6M and 
10 51 erg respectively. This leads to the final size of the sim- 
ulated PWN that is consistent with the observations. Given 



1 The full account of the convergence study will be given in the 
follow-up paper. 



the age of the Crab Nebula, it is not expected to have en- 
tered the reverberation stage and thus the interaction of the 
supernova shell with ISM is not important for the dynamics 
of the PWN. For this reason we do not introduce the ISM 
zone and continue the supernova shell solution up to the 
outer boundary of the computational domain. 

At — and — tt we impose the standard axisymmet- 
ric boundary conditions using three ghost cells. At the outer 
radial boundary, r = lOlyr, we imposed the non-reflective 
boundary conditions and at the inner radial boundary our 
boundary conditions describe the supersonic pulsar wind. 

Finally, we adopt the equation of state of ideal gas with 
the ratio of specific heats 7 = 4/3. 



2.2 Model of synchrotron emission 

The synchrotron particles are assumed to be continuously 
injected at the termination shock with the power law spec- 
trum 



f(e) = Ae- r for e < e c 



(1) 



where the cut-off energy is set to the value e c = lOOOTeV 
and the power index to F = 2.2. According to the results 
of iKennel & Coronitil (|1984b| and iVolpi et al.l 120081) this 



simplified model fits the synchrotron spectrum of the Crab 
Nebula from the optical frequencies to X-rays. In fact, the 
value of r is sensitive to the model of the nebula flow. More- 
over, the origin of radio-emitting electrons is unclear and 
the injection spectrum may have a low energy cut-off close 
to the optical energies. Although these details are important 
for the spectrum studies, they are unlikely to have strong ef- 
fect on the appearence of the nebula at any given frequency 
between optics and X-rays. The key factors here are the 
variation of magnetic field strength, and the adiabatic and 
syncrotron energy losses as they has strong effect on the 
relative brigntness of varios features. 

In order to account for the adiabatic losses (or gains) , a 
new dynamic variable, the number density of trace species n, 
is introduced in the simulations; the variation of n reflects 
the compression, or rarefaction, of the fluid element that 
carries them around. The covariant equation describing its 
evolution is 



V M (n<) = 0, 



(2) 



where is the four-velocity of the flow. In order to deduce 
the total expansion experienced by a particular fluid element 
downstream of the termination shock we also need to know 
no, the value of n when this element crossed the shock; the 
volume expansion is then given by the ratio n/no- This is 
achieved via treating no as another dynamic variable which, 
downstream of the termination shock, satisfies the transport 
equation 



V /X (n n^ /X ) = 0. 



(3) 



Upstream of the shock we assume that no equals to n and 
that n is proportional to the total energy flux, that is 



l (r ,0) = n(— ) (sin 2 6>+— ) 
\ r J ao 



(4) 



where ro, N, and ao are constants. For simplicity, we assume 
that, up to a constant coefficient, the parameter A in Eq[l] 



N 
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Figure 1. The plot shows logiop in dimensionless units. The unit 
of length in this and other plots is one light year. 



equals to no- The downstream distribution function of the 
synchrotron electrons can now be expressed as 

where Coo is the new cutoff energy due to the synchrotron 
and adiabatic losses. This is another dynamical variable 
which is evolved in the simulations according to the trans- 
port equation 



~ D 2 2 4 

-c 2 B e^n- 



(6) 



where B is the magnetic field strength as measured in 
the fluid frame and 82 — (4e 4 /9ra 4 c 7 ). In the derivation 
we assumed effective pitch-angle scattering of the electrons 
(and positrons). Upstream of the termination shock we set 
Coo = e c . The corresponding synchrotron emissivity, as mea- 
sured in the observer's frame, is then calculated using the 
approximation 



J Uoh ocn V B_ 



(?)" 



e'- r I 1 



(7) 



where v — ciB±e 2 , c\ — 3e/47rra 3 c 5 , B± is the component 
of the magnetic field normal to the line of sight in the fluid 
frame, and V = v Q h/v is the Doppler factor. 



3 RESULTS AND DISCUSSION 

Figure ^ shows the solution by the end of the simulation 
run, t = 954 yr. Compared to our previous simulations, the 
most remarkable feature of this solution is the absence of the 
equatorial disk-like outflow and the large-scale circulation in 
the nebula. Instead, we observe a rather irregular flow that 
can be described as a collection of large scale vortexes. Since 
the same behaviour is found even in the low resolution runs 
we conclude that the key factor leading to this change in the 
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Figure 2. All plots show the azimuthal component of magnetic field in CGS units, and clearly illustrate the non-linear dynamics of 
the termination shock. The left panel shows a stage at which the termination shock is fully inflated (t = 912.3 yr). The middle panel 
shows the northern shock structure being distorted (t = 926.6 yr). The right panel shows the shock complex to be highly compressed 
(t = 933.1 yr). 




properties of the PWN flow is the relaxation of the equato- 
rial symmetry condition. These vortexes pick up mass from 
the interface between the PWN and the supernova shell and 
transport it into the nebula's interior, creating the intricate 
fine structure of the density distribution. The long "fingers" 
protruding inside the nebula along the polar axis are pro- 
duced by the back- flow of the polar jets contaminated with 
the entrained mass of the supernova shell. In the proper 
three-dimensional case, we do not expect the polar jets to 
remain so straight and collimated as in our simulations. As 



the result, the backflow will be much less pronounced and 
coherent. 

Figure [2] illustrates the highly complex dynamics found 
in the central part of the simulated PWN. The structure of 
the termination shock is highly variable due to the non-linear 
interactions with the flow inside the nebula. On one hand, 
both the size and the shape of the shock change constantly 
due to "collisions" with the large amplitude waves and vor- 
texes of the nebula. On the other hand, these changes in the 
shock structure result in highly variable outflow from the 
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Figure 4. The left panel shows the fragment of time-series data for the comoving magnetic field (in arbitrary units) at the reference cell 
with coordinates r = 1.75 lyr, = 87°. Time is measured from the start of sampling. The middle panel shows the structure function of 
this time-series. The right panel shows the real part of its wavelet transform using the Morlet mother wavelet; the dashed line corresponds 
to the cone of influence. 



shock. Thus, the variations of the termination shock struc- 
ture and the variability of the flow in the nebula feed on 
each other and at this stage it is quite impossible to distin- 
guish between the cause and the effect (this is analogous to 
the famous "hen and egg" problem ) . This behaviour of the 
PW termination shock is reminiscent of that of the accre- 
tion shock in problems involving supersonic accretion (e.g. 



Blondin et al. 2003; Buras et al. 2006 Scheck et al. 2008) and 



suggests the possibility of common origin. The characteristic 
time scale of the large scale variations is around 1-2 years. 
This makes sense as the length scale of the termination shock 
is about one light year and the speed of magnetosonic waves 
is close to the speed of light. 

In contrast to the PW, the equatorial symmetry in the 
PWN is completely broken and often the fluid elements 
blown by the wind into the southern hemi-sphere (in Fig[2] 
they are coloured in blue) end up in the northern hemi- 
sphere, and the other way around. Similarly strong break- 
downs of equatorial symmetry have been observed in the 



simulations of stellar collapse (e.g. Scheck et al. 



2008 



Komis- 



sarov|2009 ) . This seems to be a general rule for the problems 



where the shocked plasma is confined to a finite or slowly 
expanding volume. 

As expected, and in agreement with previous simula- 
tions, the predominant motion near the equatorial plane is 
an outflow. The typical speed of this outflow inside the in- 
ner 2 lyr is around 0.6c. The outflow is highly inhomogeneous 
with regions of strong magnetic field following regions of rel- 
atively weak field. Such, strong variations in magnetic field 
strength lead to strong variations of synchrotron emissivity 
and ultimately to the phenomenon of expanding wisps (see 
Figure [3}. Further out, but well before reaching the super- 
nova shell the outflow slows down and mixes with the rest 
of the PWN. 

Above the equatorial plane a backflow with superim- 
posed vortexes can also be seen. This seems to be the rea- 
son for the contracting wisps occasionally observed in these 
simulations. Closer to the symmetry axis, with r < 0.5 lyr, 
another backflow is seen. Like in the previous simulations, 
this flow originates from the highly-magnetized layers of re- 



cently shocked plasma of the pulsar wind. The strong mag- 
netic hoop stress in these layers stops the outflow and turns 
it back towards the axis. The resultant axial compression 
drives the polar jets in a fashion reminiscent of the tooth- 
paste flow. The axial pinch is rather nonuniform and the 
flow structure at the jet base is highly variable. One can in- 
terpret this variability as the development of the "sausage" 
instability enhanced by the fact that the backflow is already 
highly inhomogeneous and the instability is in the non-linear 
regime from the start. The converging magnetosonic com- 
pression waves, which originate in the dynamically active 
region around the equatorial outflow, also contribute to this 
process. In the synthetic synchrotron maps this region often 
shows relatively compact and bright variable features (see 
Figure [3| which look similar to the Crab's "sprite" ([Hester 
[etaLl[2002t . 

The variability of the Crab Nebula emission was inves- 



tigated by|Hester et al. (2002) via subtraction of two images 



separated by approximately 109 days. In particular, an out- 
ward moving wisp appears on the resultant image as a dark 
wisp followed by a bright one. In Fig [3] we apply the same 
technique to our synthetic optical images. The dynamics of 
simulated wisps is remarkably similar to the that of the real 
ones. In agreement with the observations, the typical appar- 
ent speed of the wisps is c± 0.5c within the inner c± 2 lyr (In 
the mapping procedure we ignore the time-delay effect.). 
Further out the apparent speed of the wisps significantly 
decreases and we often observe wisp mergers like those de- 



scribed in Hester et al. (2002). Such behaviour agrees with 



the interpretation of wisps as magnetic inhomogeneities ad- 
vected by a nonuniform decelerating subsonic flow. 

Having observed the wisp production in our RMHD 
simulations, we then quantitatively analysed the nature of 
the nebula variability in order to establish whether the shock 
dynamics and associated synchrotron emission show any pe- 
riodicity, and if not, to determine the statistical nature of 
the fluctuations. A number of time-series were constructed 
by measuring fluid parameters at the point with coordinates 
(r,0) = (1.75 lyr, 87°) after each local time-step. The left 
panel of Fig [4] shows part of such time series for the co- 
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moving magnetic field. The middle panel of Fig|4] shows the 
structure function of the time-series. One can see that it 
has a plateau at long time lags and a broken power law at 
short time lags. The hard power law at time lags r < 0.08 yr 
is clearly a numerical noise as this time scale corresponds 
to the numerical resolution of our simulations. The softer 
power law at time lags r = 0.1 — 1 yr corresponds to well 
resolved scales and presumably captures the MHD turbu- 
lence developed in the simulations. The transition to plateau 
at r c± 1.0 yr indicates the characteristic time scale for the 
large scale process that drives the turbulence. The only can- 
didate for such a process is the oscillation of the termination 
shock. This oscillation cannot occur on scale below the light 
crossing time of the shock and this seems to explain the 
characteristic time scale revealed by the structure function. 

In order to investigate this issue further we also applied 
the method of wavelet transform. In the right panel of Fig|4] 
the real part of the wavelet transform for the Morlet mother 



wavelet (Torrence & Compo 1998) is shown. One can see 



clear evidence of quasi-periodic oscillations with the quasi- 
period P — 1.5 — 3yr. The analysis of time series for other 
flow parameters gives similar results. 

The available variability data for the Crab Nebula do 
not cover sufficiently long time period and are too scarce for 
rigorous statistical analysis. Therefore, we cannot yet make 
proper comparison between the observations and our sim- 
ulations. The observations by |Hester et al.| ([2002 ) revealed 
two wisps produced within one year. This may indicate a 
somewhat shorter characteristic timescale for wisp produc- 
tion compared to other results, but may also be just a local 
variation. Under the circumstances, we conclude that the 
MHD model explains not only the structure of the Crab 
Nebula but also its variability reasonably well. 
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